#geno<-read.table("/data1/bsi/BORA_processing/devel/eqtl/parallelize_genotypingsurvival/input_withcorrecthead.tped",sep=" ",head=T)
setwd("/data1/bsi/BORA_processing/devel/eqtl/parallelize_genotypingsurvival/1")
hazardrank<-read.table("hazardrank.txt",sep="\t",head=T)
clin<-read.delim("/data1/bsi/BORA_processing/devel/eqtl/parallelize_genotypingsurvival/clinical" ,sep="\t",quote="\"",dec=".",fill=TRUE,comment.char="", header=TRUE, col.names=c("SUBJECT","AGE","GENDER","OS","PFS","DFS","ALIVEOS","ALIVEPFS","ALIVEDFS","SITE"),as.is=c(TRUE,TRUE,FALSE,TRUE,TRUE,TRUE,TRUE,TRUE,TRUE,FALSE),colClasses=c("character","numeric","factor","numeric","numeric","numeric","integer","integer","integer","factor"))
clin$SUBJECT<-gsub("-",".",clin$SUBJECT)
tnum<-1
while(tnum<ncol(hazardrank)+1)
{
	dt <- merge(clin,t(hazardrank[tnum,]),by.y="row.names",by.x="SUBJECT")
	#dt1<-dt[sort.list(dt[,4]), ]
	dt1<-dt[rev(sort.list(dt[,11])), ]

	i<-nrow(dt1)
	#order<-order(dt1[,11])
	hazard_rank<-dt1[,11]
	N<-nrow(dt1)
	survival_time<-dt1[,4]
	k<-0
	SensArr<-c()
	oneMSpecArr<-c()
	nsens<-1
	SensArr[nsens]<-0
	oneMSpecArr[nsens]<-0
	nsens<-nsens+1
	while(i>0)
	{
		if(k != dt1[i,11])
		{
			#pi<-i
			t<-dt1[i,4]
			#coln<-colnames(geno)
			c<-hazard_rank[i]
			Pr_M_ge_c<-length(which(hazard_rank>=c))/N
			Pr_T_le_t<-length(which(survival_time<=t))/N	
			Pr_T_le_t_given_M_ge_c<-length(which(survival_time<=t & hazard_rank>=c))/N
			Pr_T_gt_t_given_M_lt_c<-length(which(survival_time>t & hazard_rank<c))/N
			Sens<-(Pr_T_le_t_given_M_ge_c*Pr_M_ge_c)/Pr_T_le_t
			if(Pr_T_le_t == 1) {
				Spec<-0
			} else {
			        Spec<-Pr_T_gt_t_given_M_lt_c*(1-Pr_M_ge_c)/(1-Pr_T_le_t)
			}
			oneMSpec<-1-Spec
			SensArr[nsens]<-Sens
			oneMSpecArr[nsens]<-oneMSpec
			nsens<-nsens+1	

			#print(c("Y",y,"dx",dx))
			k<-dt1[i,11]
		}
		i<-i-1
		
	}
	SensArr[nsens]<-1
	oneMSpecArr[nsens]<-1
	AUC<-0
	prev_sens<-0
	prev_oneMSpec<-0
	sorted_sens<-order(oneMSpecArr)
	i<-1
	while(i<=nsens) {	
		Sens=SensArr[sorted_sens[i]]
		oneMSpec=oneMSpecArr[sorted_sens[i]]
		dx<-oneMSpec-prev_oneMSpec
		y<-(Sens+prev_sens)/2
		AUC<-AUC+y*dx
		#print(c("Y",y,"dx",dx))
		prev_oneMSpec<-oneMSpec
		prev_sens<-Sens
		i<-i+1
	}
	#print(survival_time)
	#print(hazard_rank)
	#print(c("Y",y,"dx",dx))
	print(c("AUC",AUC))
	tnum =tnum+1
}
